Application of the dual-kinetic-balance sets in the relativistic many-body problem of 

atomic structure 
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The dual-kinetic-balance (DKB) finite basis set method for solving the Dirac equation for 
hydrogen-like ions [V. M. Shabaev et al., Phys. Rev. Lett. 93, 130405 (2004)] is extended to 
problems with a non-local spherically-symmetric Dirac-Hartree-Fock potential. We implement the 
DKB method using B-spline basis sets and compare its performance with the widely-employed ap- 
proach of Notre Dame (ND) group [W.R. Johnson and J. Sapirstein, Phys. Rev. Lett. 57, 1126 
(1986)]. We compare the performance of the ND and DKB methods by computing various properties 
of Cs atom: energies, hyperfine integrals, the parity-non-conserving amplitude of the 6s 1; /2 — 7si/2 
transition, and the second-order many-body correction to the removal energy of the valence elec- 
trons. We find that for a comparable size of the basis set the accuracy of both methods is similar 
for matrix elements accumulated far from the nuclear region. However, for atomic properties deter- 
mined by small distances, the DKB method outperforms the ND approach. In addition, we present 
a strategy for optimizing the size of the basis sets by choosing progressively smaller number of basis 
functions for increasingly higher partial waves. This strategy exploits suppression of contributions 
of high partial waves to typical many-body correlation corrections. 
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I. INTRODUCTION 



Applications of perturbation theory in quantum me- 
chanics require summations over a complete set of states 
of the lowest-order Hamiltonian. Usually, the relevant 
spectrum is innumerable. In practical applications such 
eigenspectra are often modeled using finite basis sets, 
chosen to be numerically complete. Since the sets are 
finite, the otherwise infinite summations become amend- 
able to numerical evaluations. 

The use of a finite basis set composed of piecewise poly- 
nomials, so-called B-splines has proven to be par- 
ticularly advantageous in atomic physics and quantum 
chemistry applications [2|. In this approach, an atom 
is placed in a sufficiently large cavity and the atomic 
wavefunctions are expanded in terms of the underlying 
B-spline set. Further, the variational Galerkin method is 
invoked and the solution of the resulting matrix eigen- 
value problem produces a quasi-spectrum for the atom. 
In non-relativistic calculations, the lowest-energy orbitals 
of the resulting basis set closely agree with those of the 
unperturbed atom, and calculations of various properties 
of the low-lying states can be carried out. In particular, 
one could generate single-particle orbitals in some suit- 
able lowest-order approximation, and use the resulting 
basis set in applications of many-body perturbation thc- 
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ory (MBPT) . 

Application of the outlined approach to the relativis- 
tic problems brings in a complication- the appearance 
in the atomic quasi-spectrum of non-physical "spurious" 
states. These states appear in the solution of the single- 
particle radial Dirac equation for k > angular symme- 
try, j =1 — 1/2 ( P1/2, ^3/2 7 ■ ■ ■ orbitals). The spurious 
states rapidly oscillate and, moreover, spoil the mapping 
of the generated quasi-spectrum onto the low-energy or- 
bitals of the atom. At the same time they are required 
for keeping the set complete. The problem of spurious 
states was discussed in the literature in details, see e.g., 
Ref. [1, |3, H, E|, and several solutions were proposed. In 
the pioneering applications of the B-splines in relativistic 
many-body problem by the Notre Dame group, Johnson 
et al. Q added an artificial potential spike centered at 
the origin to the Hamiltonian matrix. The overall effect 
was to shift the spurious states to higher energies thus 
restoring the low-energy mapping to the physical states. 
We will refer to the sets generated using this prescrip- 
tion as the Notre Dame (ND) sets. Another solution was 
to use "kinetically-balanced" sets which related the 
small and the large components of the basis set functions 
via the Pauli approximation. Recently, an extension of 
this method was proposed in Ref. [?|]- Here, due to addi- 
tional relations between the small and large components, 
the negative (Dirac sea) and positive energy spectra are 
treated in a symmetric fashion. To emphasize this built- 
in symmetry, the authors refer to their method as the 
"dual-kinetic-balance" (DKB) approach. In both meth- 
ods (by contrast to the ND prescription), the spurious 
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states were shown to be completely eliminated from the 
quasi-spectrum. 

Motivated by the success of the DKB method in com- 
puting properties of hydrogen- like ions 0,0], here we in- 
vestigate the suitability of the DKB method in model- 
ing the spectrum of the (non-local) Dirac-Hartree-Fock 
(DHF) potential. We compare the performance of the 
ND and DKB methods by computing various proper- 
ties of Cs atom: single-particle energies, hyperfine in- 
tegrals, and the parity-non-conserving amplitude of the 
6S1/2 — 7 s i/2 transition. We find that for properties in- 
volving matrix elements accumulated near the nucleus, 
the DKB method outperforms the ND approach. Other- 
wise, if the electronic integrals are accumulated far from 
the nucleus, both methods produce results of a similar 
quality. 

In addition, we investigate a possibility of using vary- 
ing number of basis set functions for different angular 
symmetries. Summations over intermediate states in ex- 
pressions of perturbation theory are carried out both over 
angular quantum numbers k and for fixed k over ra- 
dial quantum numbers. Usually, as |k| (and i) increases, 
the correlation corrections due to higher partial waves 
become progressively smaller. Intuitively one expects 
that for higher partial waves it would be sufficient to use 
smaller radial basis sets of lesser quality. This would 
reduce storage requirements for many-body calculations 
(for example, in implementing coupled-cluster formalism 
) and would speed up numerical evaluations. While such 
an approach is common in quantum chemistry see, e.g., 
Ref. [8| , the question of building the optimal B-spline ba- 
sis set was not addressed yet in relativistic many-body 
calculations. We illustrate optimizing the basis sets by 
computing the second-order energy correction for several 
states of Cs. 

This paper is organized as follows. First we recapit- 



ulate the Galerkin-type approach to generating a finite- 
basis set quasi-spectrum for the Dirac equation in Sec- 
tion [TTJ The variational method is invoked for relativistic 
action and the problem is reduced to solving the gen- 
eralized eigenvalue problem in Section III Al The DHF 
potential is specified in Section Hi B I Further we specify 
ND and DKB sets in Section III CI and boundary condi- 
tions in Section III Dl A numerical analysis of Cs atom 
is provided in Section IIIII In Sections IIII Al and IIII Bl 
we compare the performance of the ND and the DKB 
sets by computing single-particle energies, hyperfine in- 
tegrals (Section IIII A[) . and parity-nonconserving ampli- 
tudes f Section IIII B[) using both methods. The spurious 
states arising from the ND method are examined in Sec- 
tion [IlTCl In Section fill D I we consider second-order en- 
ergy corrections in the DKB method and discuss a strat- 
egy of optimizing the size of the basis set. 

II. PROBLEM SETUP 

We are interested in solving the eigenvalue equation 
Hbu (r) = e u (r) for the Dirac Hamiltonian, Ho = 
ca ■ p+[3c 2 + V nuc (r) + UDHF (r), where U nuc is the nuclear 
potential and Udhf is the mean-field (Dirac-Hartree- 
Fock) potential. Udhf is in general a non-local poten- 
tial. Assuming that both potentials are central one may 
exploit the rotational invariance to parameterize the so- 
lutions as 

,. r ,\ _ 1 f iPn K (r) Km (f) \ m 
UnK(T) ~ r { Q„ K (r) n_ Km (r) J ' [i) 

with £l Km (r) being the spherical spinors. The solutions 
depend on the radial quantum number n and the angular 
quantum number k — (l—j) (2j + 1). The large, P nK , and 
small, Q U K, radial components satisfy the conventional 
set of radial Dirac equations 



(U nuc (r) + Udhf (r) + c 2 ) P„ K (r) + c [j^--j Qn* M = e nK P nK (r) , 
-c + P nK (r) + (U n uc (r) + Udhf (r) - c 2 ) Q nK (r) = e nK Q nK (r) . 
These radial equations may be derived by seeking an extremum of the following action Q 

S = \cf {P{r)6-Q{r)-Q{r)6 + P{r)}dr+ 1 - J (P (r) , Q (r)) U DHF ( q ) * 

+ ~ J U nuc (r) (P (rf + Q (r) 2 ) dr - c 2 J Q (r) 2 dr - e i J (p (rf + Q (r) 2 ) dr + AS bnd + AS s ^ r , (2) 

I 



where the K-dependent Pauli operators are defined as 

dr r 



Radial integrals here and below have implicit limits from 
r = to r = R, where R is the radius of the confining cav- 
ity. Boundary conditions may be imposed by adding the 
term AS' bnd to the action. The term AS' spur controls the 
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appearance of the spurious states in the quasi-spectrum. 
We will specify these two terms below. 



A. Reduction to the matrix form 

We employ two finite basis sets {k (r)} and {si (r)} , 
i = 1,27V, which, since we are interested in solving the 
(angularly-decoupled) radial equations, may depend on 
the angular quantum number k. We may expand the 
large and small components in terms of these bases 



2 A' 



P(r) = ^2 Pi k(r), 

i=l 
2N 



(3) 



the expansion coefficients being the same for both the 
large and the small components. The above expansions 
are plugged into the action, Eq. and, further, its ex- 
ternum is sought by varying the expansion coefficients. 
As a result one arrives at the following generalized eigen- 
value equation 



Ap = eBp, 



(4) 



where A and B are 2N x 2N matrixes and p is the vector 
of expansion coefficients in Eq.([3]). The matrix elements 
A are given by 



Aij — Du + Vij 2c Sij + A 



bnd _,_ ^spur 



(5) 



The matrixes entering the definition of A correspond to 
various pieces of the radial Dirac equations, 

Aj = c (| h (r) ± Sj (r) dr-Jk (r) (^) 8j (r) dr) 

+ c (J h (0 ^ ( r ) dr - J h ( r ) (~) s i ( r ) dr ^j > 

(6) 



Vij = / Vnuc (r) [U (r) lj (r) + s z (r) Sj (r)] dr 



DHF ) 



Sij = Si (r) Sj (r) dr, 



(7) 
(8) 



with the matrix elements of the DHF potential, 

(VbHp)ij, S iven below - The terms A u" d and 
arise from the boundary and "spurious state" correc- 
tions in the action. Finally, the matrix B, Bij = 
J [h ( r ) h ( r ) + s i ( r ) s j ( r )] dr, reflects the fact that the 
basis sets may be non-orthonormal. 



B. Potentials 

The nuclear potential V nuc (r) is generated for a nu- 
cleus of a finite size. We employ the Fermi distribution 



with the nuclear parameters taken from Ref. [9]. As to 
the Dirac-Hartree-Fock potential, we employ the frozen- 
core approximation. In this method, the calculation is 
carried out in two stages. First, the core orbitals are 
computed self-consistently. Second, based on the pre- 
computed core orbitals, the DHF potential is assembled 
for the valence orbitals and the valence orbitals are de- 
termined. In the valence part of the problem, the core 
orbitals are no longer adjusted. Explicitly, for a set of 
the angular symmetry k, 



h Si ) Vdhf 



= (^e>hf) y + (VgpSOy , (9) 

aEcore 

x J v (a, a, r) [k (r) lj (r) + s l (r) Sj (r)] dr, 

(10) 



(Wf), 



aGcore L 

x J v L (a, j, r) [U (r) P a (r) + s t (r) Q a (r)] dr, 

(11) 

with the conventionally defined multipolar contributions 



v L (b, a, r) = / -£r [P a (/) P fc (/) + Q a (/) Q b (r')] dr' , 



and the angular coefficient 



A 



K a LKu 



ja jb L 

-1/2 1/2 



(12) 



(13) 



C. ND and DKB sets 

Now we specify the Notre Dame (ND) and the dual- 
kinetic-balance (DKB) basis sets. Both operate in terms 
of B-spline functions and first we recapitulate the relevant 
properties of these splines. A set of n B-splines of order k 
is defined on a supporting grid {U} , i = 1, n + k. Usually, 
the gridpoints are chosen as 



f 11 



ta = 
t n +i 



tk = 0, 

— Ui+k = R- 



In our calculations the intermediate gridpoints are dis- 
tributed exponentially. B-spline number i of order k , 
B^ 1 (r), is a piecewise polynomial of degree k — 1 inside 
U < r < tj+fc. It vanishes outside this interval. This 
property simplifies the evaluation of matrix elements be- 
tween the functions of the basis set. In addition, we will 
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make use of the fact that as r — > 0, the first k splines 
behave as (all the remaining splines are zero) 



(14) 



and at r = R, all splines vanish except for the last spline, 



B 



(k) 

i—n ' 

The Notre Dame set is defined as 

/ Pf ) (r) 1 < i < n 



n < i < 2n ' 
1 < i < n 
B ( ^ n (r) n<i<2n ' 



(15) 



It corresponds to an independent expansion of the large 
and small radial components into the B-spline basis. The 
DKB set involves the Pauli operators and enforces a "ki- 
netic balance" between contributions to the components: 




Bf (r) 



Ki<n 



_ / hP% B i k) ( r ) ^<i<n 



(16) 



B\% (r) n < i < 2n 

Notice that, as discussed below, to satisfy the boundary 
conditions, we will use only a subset of the entire DKB 
basis. 



D. Spurious states and boundary conditions 

With the ND set, the spurious states are shifted 
away to the high-energy end of the quasi-spectrum by 
adding the following action @ (A5* spur ) ND = f P (0) 2 - 
§ P (0) Q (0) for k < and (AS s P ur ) ND = c 2 P (0) 2 - 
|P (0) Q (0) for k > 0. This correction may be seen as 
arising from an artificial 5— function-like potential cen- 
tered at the origin. Unfortunately, as shown below in 
numerical examples, this additional spike perturbs the 
behavior of the orbitals near the nucleus. (As discussed 
below (AS spur ) ND also sets the boundary conditions at 
r = 0.) The DKB set does not have the spurious states 

at all, so that (AS ,spur ) DK B = °- 

We need to specify boundary conditions at r — and 
at the cavity radius, r = R. We start by discussing the 
boundary conditions at r — 0. For a finite-size nucleus 
the radial components behave as 



P nK oc r l+1 and Q nR oc r l+2 for k < , 
P nK oc r l+1 and Q nK oc r l for n > . 



(17) 



In the Notre Dame approach, the boundary conditions 
are imposed variationally by adding the boundary terms 
to the action integral. Varying AS' spur , Ref.fjf, effec- 
tively reduces to P (0) = 0. In practice, because of the 
variational nature of the ND constraint, the large com- 
ponent, although being small, does not vanish at the 



origin, and the limits, (|17p. are not satisfied. Alterna- 
tively, Froese-Fischer and Parpia [l(| proposed to im- 
pose P (0) = Q(0) = by discarding the first B-splinc 
of the set (this is the only B-spline that does not van- 
ish at r = 0). This is a "hard" constraint, since the 
wavefunction, Eq.([3]), would vanish identically at the ori- 
gin. In our calculations, because of our motivation in 
building the smallest possible basis set, we extend this 
scheme further. We exploit the power-law behavior of 
the B-splines, Eq.(fT4|). and match it to the small-r limits 
(fiTf . To satisfy the matching, we need to include the 
B-splines starting with the sequential number (i m i n must 
be smaller than the order of the splines k). 



For S1/2 (« = -1) and p 1/2 (« = +1) 



i n = 2, and this 

is equivalent to the boundary condition of Ref. [Ioj |. 
For higher partial waves, however, an increasingly 
larger number of splines is discarded: e.g., for 
/V/2 ( K — — 4) , i m ; n = 5. One should notice that for a 
basis that includes partial waves t < £ max , for a faithful 
representation of the small-r behavior in all the partial 
waves, one needs to require the order of the splines to 
be at least of k = £ max + 3. In particular, for k = 7, 

^max — 4. 

When the first B-spline of the set, P} =1 (r), is included 
in the basis (as in the ND approach), there is another dif- 
ficulty in the calculations: since it's value does not vanish 
at r = 0, the matrix elements Di^ n +i and D^i, which 
contain matrix elements of 1/r, arc infinite in absolute 
value [6]. In practical calculations, one uses Gaussian 
quadratures to evaluate this integral, so the result of the 
integration is finite. Yet this introduces arbitrariness in 
the ND calculations and may be a reason for a relatively 
poor representation of the orbitals near the origin. 

At the cavity radius, to avoid overspecifying the 
boundary conditions for the Dirac equation, the ND 
group used the boundary condition P (R) = Q (R). As 
with the conditions at the origin, this relation was "en- 
couraged" variationally. In our calculations (similarly to 
Ref. 0) we use the "hard" condition P (R) = Q (R) = 
by removing the last B-spline from the set. Examination 
of the resulting orbitals reveals that the wavefunctions 
acquire a non-physical inflection towards the end of the 
supporting grid, while the ND orbitals behave properly. 
Further numerical experimentation, however, shows that 
the inflection does not degrade the numerical quality at 
least for the atomic properties of the low-lying bound 
states of interest. At the same time, throwing away the 
last B-spline reduces the number of basis functions and 
leads to a more compact set. 

To summarize, we will use the DKB basis set that 
includes B-splines with sequential numbers « m in('*)) 
Eq. (|18p to i max = n — 1. We will simply refer to this 
choice as the DKB basis. When we refer to N = 40 
DKB functions for a given partial wave, it would imply 
larger underlying B-spline set, e.g., for s-waves the total 
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number of B-splines would be n — 42. For the ND basis 
11 = N, as all B-splines participate in the expansion. 

In the remainder of this paper we present results of 
numerical analysis for 133 Cs atom. It is an atom with 
a single valence electron outside a closed-shell core. As 
the first step we carry out finite-difference Dirac-Hartree- 
Fock calculations for Cs core. The core orbitals are fed 
into the spline code where they are used to compute the 
matrix elements of the Vdhf potential, Eq.©. As in 
Ref. [B[ the numerical accuracy is monitored by compar- 
ing the resulting quasi-spectrum with the DHF energies 
from the finite-difference code. 



III. NUMERICAL EXAMPLES FOR CS ATOM 

Here we provide numerical examples involving both 
ND and DKB sets for Cs atom. In the two Sections im- 
mediately following we compare the performance of the 
ND and the DKB sets. We generate the quasi-spectrum 
using both ND and DKB sets and carry out compar- 
isons for single-particle energies and hyperfine integrals 
in Section IIII Al and parity- nonconserving amplitude in 
Section IlII Bl Section [ill CI contains an analysis of spuri- 
ous states in the ND method. In Section llll Dl wc analyze 
second-order energy corrections and discuss a strategy of 
optimizing the size of the basis. 



A. Energies and hyperfine integrals 

We compare ND and DKB quasi-spectrums with ener- 
gies obtained using a finite-difference DHF code for the 
low-lying valence states in Table [J] The ND and DKB 
calculations were carried out using N — 40 basis func- 
tions for B-splines of order k — 7. We used a cavity of 
radius R = 50 bohr. For the cavity of this size, only a few 
lowest-energy orbitals remain relatively unperturbed by 
the cavity. From examining Table HI it is clear that both 
ND and DKB sets have a similar accuracy for energies. 

In the second part of Table U we compare values of the 
radial integrals entering matrix elements of the hyperfine 
interaction due to the electric (EJ) and magnetic (MJ) 
multipolar moments of a point-like nucleus 



Iej (n«) 



dr 



T {PL (r) + Ql K (r)) 



dr 



I M j (n/s) = 2 / -^iPnK (r) Q nK (r) . 

The angular selection rules require j < J/2. We use iden- 
tical integration grid for all three cases (DHF, ND, DKB) 
listed in the Table. The grid is sufficiently dense near 
the origin, it contains about 100 points inside the nu- 
cleus. The numerical integration excludes the first in- 
terval of the grid. From examining the Table we find 
that the DKB set outperforms the ND basis. While 
the ND set still recovers two-three significant figures for 



the magnetic-dipole coupling, it produces wrong results 
for electric-quadrupole and magnetic-octupole integrals. 
Certainly, the accuracy in the ND case improves for a 
larger basis set, but larger basis sets come at an addi- 
tional computational cost. We carried out similar com- 
parisons for matrix elements of the electric-dipole oper- 
ator. As for the energies, we find that both the ND and 
DKB sets perform with a similar numerical accuracy. 

As we have mentioned, (AS' spur ) ND variationally en- 
courages the boundary condition P (0) = 0. However, 
there is no such explicit encouragement for Q(0). Here 
we qualitatively discuss the observed properties of the ra- 
dial components near the origin for the Cs ND set. We see 
that, though their general behavior is to approach zero, 
the large component wavefunctions often have small im- 
proper inflections or oscillations near the origin. Small 
component wavefunctions, on the other hand, often do 
not even approach zero. These improper behaviors seem 
to be exemplified as we look at states higher in the spec- 
trum. As we have seen here, such improper behavior of 
both the large and small component radial functions near 
the origin prove detrimental for properties accumulated 
near the nucleus. We conclude that while producing the 
results of a similar quality for matrix elements accumu- 
lated far from r = 0, the DKB set provides a better 
approximation to the atomic orbitals near the nucleus. 



B. Parity- nonconserving amplitude 

So far we examined properties of the individual ba- 
sis set orbitals with sufficiently low energies. The real 
power of the finite basis set technique lies in approxi- 
mating the entire innumerable spectrum by a finite size 
quasi-spectrum. This is important, for example, in com- 
puting sums over intermediate states (Green's functions) 
in perturbation theory. 

From the discussion of the preceding Section it is clear 
that the difference in quality of the ND and DKB basis 
sets is expected to become most apparent for the prop- 
erties accumulated near the nucleus. Here, as an illus- 
trative example, we consider the parity-nonconserving 
(PNC) amplitude for the 6S1/2 — > 7<Si/2 transition in 
133 Cs. This amplitude appears in the second order of 
perturbation theory for the otherwise forbidden dipolc 
transition and it can be represented as a sum over inter- 
mediate states rap 1/2 



(7si /2 \D\npi/2)(npi/ 2 \H w \6si/ 2 ) 



£6s 1/2 e np 1 



-Epnc — 

n=2 

+ ( 7s i/^ H w\ n Pi/2){npi/2\D\6s 1/2 ) ^ 

n=2 £ 6s i/2 — £ np 1/2 

Here D and i?w are electric-dipole and weak interaction 
(pseudo-scalar) operators, and Ei are atomic energy lev- 
els. We will compute this expression in the single-particle 
approximation. Specifically, the index ra runs over the en- 
tire DHF quasi-spectrum for py 2 partial wave, including 
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State Set 


Energy 




Ml HFI 




E2 HFI 


M3 HFI 


6si/ 2 FD 


-0.1273680 




1.114751[- 


-11 

J 






DKB 


-0.1273674 




1.1 14741 [- 


-1] 






ND 


-0.1273682 




1.121812[- 


-1] 






7s 1/2 FD 


-0.5518735[- 


-11 

J 


3.063077[- 


-21 

J 






DKB 


-0.5518714[- 


-1] 


3.063069[- 


-2] 






ND 


-0.5518750[- 


-1] 


3.084164[- 


-2] 






6pi 10 FD 


-0.8561589[- 


-H 


-1.252026[- 


-21 






DKB 


-0.8561576[- 


-1] 


-1.252018[- 


-2] 






ND 


-0.8561616[- 


-1] 


-1.218362[- 


-2] 






6D3/2 FD 


-0.8378548[- 


-11 

J 


4.649107[- 


-31 

■ J 


6.693978[-l] 


8.725496(0] 


DKB 


-0.8378543[- 


-1] 


4.649117[- 


-3] 


6.694037[-l] 


8.744759(0] 


ND 


-0.8378538[- 


-1] 


4.649454[- 


-3] 


1.591405[+2] 


1.924786(7] 


5d*/n FD 


-0.6441964[- 


-11 


-3.543808[- 


-31 

- J 


1.702467[-1] 


-8.950592(1] 


DKB 


-0.6441961[- 


-1] 


-3.543799[- 


-3] 


1.702736[-1] 


-9.202786(1] 


ND 


-0.6441970[- 


-1] 


-3.038702[- 


-3] 


1.643100[+5] 


3.663198(10] 


5d 5/2 FD 


-0.6452977[- 


-1] 


2.257618[- 


-3] 


1.562954[-1] 


1.938370(1] 


DKB 


-0.6452976[- 


-1] 


2.257616[- 


-3] 


1.562953[-1] 


1.938705(1] 


ND 


-0.6452969[- 


-1] 


2.257607[- 


-3] 


7.144243[+0] 


7.235873(5] 



TABLE I: Comparison of the energies and radial integrals of the hyperfine interaction in the Dirac-Hartree-Fock approximation 
for Cs. FD marks values produced by a finite-difference code. DKB and ND values are generated with dual-kinetic-balance 
and Notre Dame B-spline basis sets. In both cases we used N = 40 basis functions for B-splines of order k = 7 in a cavity of 
R — 50 bohr. Notation x[y] stands for x x 10 y . 

Further insights may be gained from examining in- 
dividual contributions of the intermediate states in the 
PNC amplitude. We plot individual contributions of the 
intermediate state np\/2 to the PNC amplitude in Fig. [I] 
Both terms in Eq. (fT9")) are included. We computed the 
data using N = 100, k = 11 ND and DKB basis sets 
generated in a cavity of R = 50 bohr. From the plot, 
we observe that the dominant contribution arises from 
the low- lying valence states. As n increases, the contri- 
butions become quickly suppressed (there is 10 orders of 
magnitude suppression for n w 50). This is due to both 
increased energy denominators and decreased density at 
the nucleus for high n. Comparison between the basis 
sets reveals that their contributions are identical until 
n = 17. For higher principal quantum numbers, the DKB 
contributions monotonically decrease, while ND contri- 
butions become irregular. Moreover, at n = 28 the ND 
contributions start to flip signs with increasing n. Gener- 
ally, this oscillating pattern would lead to a deterioration 
of the numerical accuracy. We believe that the described 
irregularity is again due to the aforementioned improper 
behavior of orbitals near the nucleus, as the matrix ele- 
ments of the weak Hamiltonian is accumulated largely in 
this regime. 

To summarize, the DKB set is numerically complete 
and is well suited for carrying out practical summations 
over intermediate states in perturbation theory. The 
comparison with the ND set shows that, at least for the 
PNC amplitude, the ND convergence pattern becomes 



core orbitals. The weak Hamiltonian reads 
G F 

H w = ^|Qw/W(r)75 , (20) 

where Gf is the Fermi constant, Qw is the weak charge, 
75 is the Dirac matrix (it mixes large and small compo- 
nents), and Pnucif) is the neutron density distribution. 
For consistency with the previous calculations the p nue (r) 
is taken to be the proton Fermi distribution of Ref. [ll[ • 
Notice that the matrix elements of the weak interaction 
are accumulated entirely inside the nucleus. 

We evaluate the sum (fl"9")) using the DKB and the ND 
sets with N — 40 basis functions of order k = 7. The 
integration grids are the same in both cases and include 
large number of points (~ 100) inside the nucleus. The 
PNC amplitude is conventionally expressed in units of 
10~ n i|e|ao(— Qw/N n ), where A^ n = 78 is the number 
of neutrons in the nucleus of 133 Cs. In these units, the 
results are 

E*» c = -0.740, 

^pnc = -0.7395(A^ = 40,fc = 7), 

-^pnc = -0.8546(A^ = 40,fc = 7). 

The finite-difference value is taken from Ref. [llj . Again 
we note that the DKB set offers an improved performance 
over the Notre Dame set. Reaching the comparable ac- 
curacy in the ND approximation requires a larger basis 
set. For example, N = 75, k — 9 ND set reproduces the 
DKB result for the PNC amplitude. 
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affected by the inaccurate representation of the orbitals 
near the nucleus (this is exemplified for states higher in 
the spectrum), while the DKB set exhibits a monotonic 
convergence. Additionally, the DKB set is devoid of spu- 
rious states, and therefore the incidental inclusion of spu- 
rious states in summations over intermediate states is of 
no concern. 
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FIG. 1: (Color online) Individual contributions of the pi/2 in- 
termediate states to the PNC amplitude, Eq. (|19[) . as a func- 
tion of the principal quantum number n. The units of the 
PNC amplitude are 10 _11 i|e|ao( — Qw/N n ). The results from 
the DKB basis are represented by squares and those from the 
ND basis by circles. Both sets contain N = 100 basis orbitals 
and use identical integration grids. Due to the employed loga- 
rithmic scale, we plot the absolute values of the contributions. 
We use hollow symbols to indicate negative contributions and 
filled symbols for marking positive contributions. 



C. Analysis of spurious states in the ND method 

Here we analyze the effect that the additional term 
(AS^^^p has on spurious states which occur in the 
ND method. We start by taking 



A5 spur = 



xfP(O) 2 
xc 2 P (O f 



§P(0)Q(0), 
§P(0)Q(0), 



K < 
K > 



(21) 



with x an adjustable parameter. For the case of x = 1 
this is equivalent to (AS' spur ) ND . We mentioned pre- 
viously that (AS rspur ) ND variationally encourages the 
boundary condition P (0) = 0; this is also true for any 
choice of x. The specific choice of x effectively alters the 
degree of such variational "encouragement." 

We find that setting x = results in a single spuri- 
ous state which appears as the lowest energy eigenstate 
for each k > 0. By subsequently increasing the value 
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FIG. 2: Location of the spurious state in the spectrum of Cs 
for 1 < k < 4 at various values of x. Here s corresponds to 
the location of the spurious state in the spectrum (e.g., s = 2 
corresponds to the spurious state appearing as the second 
lowest energy state for that k). We used a set with N = 40 
B-splines of order k = 7 confined to a cavity of radius R = 50 
bohr. 



of x towards x = 1 we may then deduce the effect that 
(AS' spur ) ND has on these spurious states as well as the 
rest of the spectrum. We find that small increases in x 
from x = causes the energy of the spurious state to 
increase, while the other states remain essentially unaf- 
fected (except for the case of near degeneracy with one 
of these states; this situation is discussed below). As x is 
increased from zero, we may watch as spurious states for 
each k > first appear as the lowest energy state, then 
move up to the second lowest energy state, then to the 
third lowest energy state, etcetera. Fig. [2] displays this 
effect for the case of Cs. 

It is also interesting to analyze the effect of the spurious 
state when its energy is nearly degenerate with another 
state in the spectrum (we will refer to this other state 
as the "genuine" state). As the energies approach degen- 
eracy by varying x, we observe that the spurious state 
begins to mix in with the genuine state. The first evi- 
dence for this is that the energy of the genuine state starts 
to become affected by the presence of the spurious state. 
Secondly, we see that the radial wavefunctions P(r) and 
Q(r) of the genuine state begin to acquire non- physical 
"bumps" that oscillate in a way that corresponds to the 
rapid oscillations of the respective spurious state radial 
functions. As the energy becomes nearly degenerate, the 
two states mix to such a degree that it is not even possi- 
ble to define one as the "genuine state" and the other as 
the "spurious state." This effect is shown in Fig.[3l where 
the genuine state is taken as the 4d 3 / 2 state of Cs. As 
x is increased such that the spurious state becomes em- 
bedded in the quasi-continuum part of the spectrum, it 
mixes with multiple neighboring states, and at this point 
it becomes difficult to track or even define the spurious 
state. 

Presumably increasing x to x = 1 (the Notre Dame 
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(4d 3/2 , spur), e = -3.503 
(4d yy spur),e = -3.471 




-1.5 




FIG. 3: (Color online) Behavior of the large component radial 
wavefunction P(r) of the 4d 3 / 2 (k = 2) state of Cs when the 
spurious state is (a) just below it in the spectrum (x = 0.052), 
(b) nearly degenerate with it (a; = 0.054), and (c) just above 
it in the spectrum (x = 0.056). For the case (b), the notation 
(4^3/2, spur) is used to indicate that the states are some linear 
combination of the two states and not easily defined as one or 
the other. Here e refers to the energy in a.u. with rest energy 
subtracted off. For reference, the ND set (x = 1) gives the 
lowest three d 3 / 2 states to have energies: 3d, 



4(/ 



3/2, 



-28.310; 



3/2, 



e = -3.486; 5d 3/2 , e = -0.064. We used a set with 
N = 40 B-splines of order k — 7 confined to a cavity of 
R — 50 bohr; note that these plots only extend to r — 3 bohr, 
however. 



case) shifts the spurious states all the way to the end of 
the spectrum. Evidence for this is seen, for example, in 
the d 3/2 set used for Fig. [3] (N = 40, k = 7, R = 50 
bohr). Here the last few levels (excluding the very last 
level) have an energy spacing on the order of 10 6 a.u., 
whereas the energy spacing to the final level is on the 
order of 10 9 a.u. Furthermore, increasing x past x = 1 
only has a substantial effect on the energy of the final 



level (for x = 10, 000 the energy spacing to the final level 
is then on the order of 10 13 a.u.). 

Up to this point we have been exclusively considering 
cases with k > 0, as these have been the only angu- 
lar symmetries for which spurious states have previously 
been known to occur. Now we shall consider the effect 
that (AS' spur ) ND has on the cases of k < 0. As with the 
k > cases, we see that increasing x past x — 1 only 
has a substantial effect on the energy of the final level. 
This is an indication that a spurious state may actually 
lie at the end of the spectrum for k < angular symme- 
tries as well. In fact, we find that setting a; to a negative 
value significantly below x — results in a single spuri- 
ous state which appears as the lowest energy eigenstate 
for each k < 0. By subsequently increasing x from this 
point, we may watch as each spurious state is shifted to- 
wards the higher energy end of the spectrum, similar to 
what is observed with k > spurious states. 

Now we return briefly to subject of the matrix elements 
£>i, n +i and D n+ \,i of the ND basis, which are suspected 
to contribute to poor representation of orbitals near the 
nucleus. From Eqs. ^ and (|15p we see that these include 

the integral J ~ (r) dr, which is infinite due to 

the non-vanishing property of the first B-spline at r — 
0. Numerically this integral is evaluated by Gaussian 
quadrature, producing finite values. We observe that for 
the (N — 40, k = 7, R = 50 bohr) Cs set, increasing the 
numerical value of this integral to 20 times its Gaussian 
quadrature value results in the reappearance of spurious 
states as the lowest energy eigenstates for each n > 0. 
Simultaneously, the energy of highest energy eigenstate 
for each k < is increased substantially. Evidently the 
capability of (AS' spur ) ND to shift the spurious state to 
the end of the spectrum for k > angular symmetries is 
reliant upon the inaccurate numerical evaluation of this 
(theoretically) infinite-value integral. 

The claim made in this section of observing spurious 
states for n < angular symmetries may seem surprising 
at first. Shabaev et al. [7| have proved that an indepen- 
dent expansion of large and small radial components with 
a finite set of basis functions leads to a single spurious 
state for k > angular symmetries. This proof assumes 
the basis functions to vanish at the origin, and the result 
of this proof is consistent with experience when such ba- 
sis functions are employed. Arguably, connection is lost 
immediately with this proof because the ND set includes 
the first B-spline, which does not vanish at the origin. 
As we have seen here, the ND method depends on nu- 
merical inaccuracies in evaluating infinite- value integrals 
in order to manipulate spurious states arising in k > 
cases. Because the ND method also includes the first B- 
splinc for k < angular symmetries (and hence similar 
numerical inaccuracies) , we would have no reason to dis- 
count the possibility of spurious states from occurring in 
these cases as well. 
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D. Optimizing the basis set: second-order energy 
correction 

In the preceding sections we established that the DKB 
sets are more robust than the ND bases. For a compa- 
rable size of the basis set the DKB basis exhibits better 
numerical accuracy for properties accumulated at small 
radii. Or, we may say that for a fixed numerical accuracy, 
the DKB basis may contain smaller number of basis func- 
tions. In this section we investigate a related question: 
what the smallest possible basis set is for a given numeri- 
cal accuracy. Keeping the set as small as possible speeds 
up many-body calculations that usually require multi- 
ple summations over intermediate states. Also, smaller 
basis sets reduce storage requirements for expansion co- 
efficients in all-order techniques such as configuration in- 
teraction or coupled-cluster methods. 

Quantifying the numerical accuracy requires choosing 
some metric, which characterizes deviation of the selected 
property for a given basis from its exact value. Appar- 
ently, one should select the "metric" so that it can be 
easily computed and is related to the relevant atomic 
properties. As an example of an optimization measure, 
here we choose the second-order correction to the energy 
of a valence electron, . 

(2) 

In the frozen-core DHF approximation Ei ' is the 
leading many-body correction to the energy. It is 
given in terms of the Coulomb integrals gijki — 
J G?ld2uj(l)«j(2) (I/V12) Ufc(lVuj(2) and single particle 
energies £j as (see, e.g., Ref. [l2[ ) , 

t^(2) \ Qabvndvnab \ Svamnfjmnva 

K = — t — z — z — — z — z — z — ' 

abn mna 

(22) 

where g%jki — 9ijH ~ 9ijlk is the antisymmctrized 
Coulomb integral. The summations are carried out over 
core orbitals (labels a and b) and virtual (non-core) or- 
bitals (labels m and n). Each summation implies sum- 
ming over principal quantum numbers, angular numbers 
k, and magnetic quantum numbers. The summation over 
magnetic quantum numbers may be carried out analyti- 
cally and we are left with summations over radial func- 
tions. 

Wc define a contribution of an individual partial wave 
£, 6Ei 2) (£), as the difference SE {2 \i) = \i)-E { f' \t- 
1), where E^\i) stands for truncated Eq. (f2"2")l : it in- 
cludes summations over intermediate states (both core 
and virtual) with the orbital angular momentum up to t. 
Since the calculations necessarily involve Coulomb inte- 
grals between orbitals of different angular momenta, the 

numerical error in 8E\r'{() is affected by the accuracy of 
representation of all partial waves up to I. 

First, in Table HT1 we present results for a large set of 
N = 100 basis functions for each partial wave. We use a 
sufficiently large cavity of R = 50 bohr in this calculation. 
These results will serve as a benchmark for comparisons 
with the less complete (optimized) sets. We observe that 



the dominant (60%) contribution comes from d-waves, 
1 = 2. Qualitatively, the second-order energy correction 
may be described as core-polarization effect. The outer 
3d-orbitals of the core are relatively "soft" and are easily 
polarizable. Since the core does not contain / and higher 
partial waves, after peaking at I = 2, the partial-wave 
contributions become suppressed as I increases. Qualita- 
tively, this suppression arises due to increased centrifugal 
repulsion of higher partial waves and associated smaller 
overlaps with core orbitals in the Coulomb integrals of 
Eq. I|22p. For example, I = 6 contributes only 0.6% of 
the total value. 

Now we would like to minimize the size of the set 
by choosing a different number of radial basis functions 
N K for different angular symmetries k. To preserve 
a numerical balance between the fine-structure compo- 
nents (for example, this may be important while recov- 
ering non-rclativistic limit) we keep the same number of 
functions for a given orbital angular momentum £, e.g., 
N Pl/2 =N P3/2 . 

In Table Hit we present an example of an optimized ba- 
sis (marked as "Small set"). We also list resulting numer- 

(2) 

ical errors for each partial wave by comparing 8E &s '(l) 
with the result from the "Large set" calculations. We see 
that while the basis-set error in higher partial waves is 
as large as 40%, this hardly makes any influence on the 
total value of the correlation energy, because of contri- 
butions of higher I are suppressed. The total value of 
the correlation correction differs by about 1% from it's 
saturated value. Considering that the correlation contri- 
bution to the energy is about 10% for the 6s state, the 
less complete set would introduce only 0.1% error for the 
total ionization energy. 

So far we discussed the correlation energy correction 
for the 6S1/2 state. The optimized set remains sufficiently 
robust for other low-lying states as well. We have car- 
ried out a comparison similar to Table HT1 for 6|>i/2, 6p 3 / 2 , 
5c?3/2, and 5(^/2 states. In all these cases the difference 

(2) 

between the total Ei, values computed with the "Large" 
and "Small" sets is about 0.5%. In practical calculations 
one is often required to reproduce a number of proper- 
ties with the same set. The atomic properties may be 
quite dissimilar- like hyperfine-structure interactions ac- 
cumulated near the nucleus and dipole matrix elements 
determined by the valence region. Apparently, one has to 
carry out a similar low-order MBPT analysis for the rel- 
evant quantities to verify the suitability of the optimized 
set. 

From the computational point of view, using the op- 
timized sets speeds up the numerical evaluations. In 
our illustrative example, the "Large set" contains N t — 
(1 + 2 x 6) x 100 = 1300 orbitals, while the "Small set" 
is about four times smaller (N t = 285 orbitals). The re- 
sulting speed-up is sizable: our computations of Eq. ([22]) 
for the 6S1/2 state (I = 6) are about 14 times faster 
with the optimized set, as expected due to _/V t 2 scaling 
of the number of contribution in the most computation- 
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ally demanding second term of Eq. (j22f . Similar scaling 
should hold for storage of expansion coefficients in all- 
order methods, e.g., for storing triple excitations (l3l | one 
expects -ZV t 3 scaling of the storage size. Usually higher 
partial waves produce larger number of angular channels 
and the scaling should be even steeper than iV t 3 . Fur- 
ther speed-up in MBPT summations and reduction in 
storage size may be attained by skipping a few basis set 
functions at the upper end of the quasispectrum. Such 
additional truncation of the spectrum becomes apparent 
from Fig. [TJ where contributions to the PNC amplitude 
for the upper three-quarters (n > 25 out of N = 100) 
of the quasi-spectrum affect the total value below 10~ 4 
level of accuracy. 



1 


Larj 


?e set Small set 


Error 




(N,k) 


5Ef a \i) (N,k) 









(100,11) 


-0.0000130 (35,7) 


-0.0000122 


6% 


1 


(100,11) 


-0.0020027 (35,7) 


-0.0019936 


0.5% 


2 


(100,11) 


-0.0105623 (30,5) 


-0.0105373 


0.2% 


3 


(100,11) 


-0.0039347 (25,4) 


-0.0039095 


0.6% 


4 


(100,11) 


-0.0007563 (15,4) 


-0.0007269 


4% 


5 


(100,11) 


-0.0002737 (10,4) 


-0.0002272 


20% 


6 


(100,11) 


-0.0001182 (10,4) 


-0.0000844 


40% 


Total 




-0.0176609 


-0.0174912 


1% 



TABLE II: Contribution of individual partial waves I to the 
second-order energy correction for the ground 65*1/2 state of 
Cs atom. We use the DKB basis set of N basis functions con- 
structed from a subset of B-splines of order k (label (iV, k)). 
Calculations are carried out in a cavity of R = 50 bohr for 
two basis sets "Large" and "Small". Numerical integration 
grids are identical for both sets. The column marked "Error" 
refers to a relative error in a given partial-wave contribution, 
SE^J (£), caused by switching from the "Large" to the "Small" 
basis set. 



between hyperfine levels, and many-body corrections to 
hyperfine interactions require accurate representation of 
atomic orbitals at both small and intermediate electron- 
nucleus distances. Solving the many-body problem in 
high orders of perturbation theory additionally calls for 
an efficient representation in terms of the basis sets. The 
dual-kinetic-basis set is shown here to adequately meet 
both these demands. 

Previously the DKB basis was applied to systems with 
a single electron in a Coulomb held, i.e., hydrogen- like 
ions[a,0|. Here we extended the DKB method to many- 
electron systems by generating the single-particle quasi- 
spectrum of the Dirac-Hartree-Fock potential. Several 
numerical example for Cs atom were presented. We 
showed that the DKB method outperforms the widely- 
employed Notre Dame B-spline method Johnson et al. Q 
for problems involving matrix elements accumulated at 
small distances. 

In addition, we presented a strategy for optimizing the 
size of the basis sets by choosing progressively smaller 
number of basis functions for increasingly higher partial 
waves. This strategy exploits suppression of contribu- 
tions of high partial waves to typical many-body correla- 
tion corrections. 
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